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ABSTRACT 

Planetary Nebulae (PNe) are amongst the most spectacular objects produced by stellar evolution, 
but the exact identity of their progenitors has never been established for a large and homogeneous 
sample. We investigate the relationship between PNe and their stellar progenitors in the Large 
Magellanic Cloud (LMC) by means of a statistical comparison between a highly complete spectroscopic 
catalog of PNe and the spatially resolved age distribution of the underlying stellar populations. We 
find that most PN progenitors in the LMC have main-sequence lifetimes in a narrow range between 5 
and 8 Gyr, which corresponds to masses between 1.2 and 1.0 M 0 , and produce PNe that are visible 
for 27 ± 6 kyr. We tentatively detect a second population of PN progenitors, with main-sequence 
lifetimes between 35 and 800 Myr, masses between 8.2 and 2.1 M 0 , and average PN lifetimes of 
lllg kyr. These two distinct and disjoint populations strongly suggest the existence of at least two 
physically distinct formation channels for PNe. Our determination of PN lifetimes and progenitor 
masses has implications for the understanding of PNe in the context of stellar evolution models, and 
for the role that rotation, magnetic fields, and binarity can play in the shaping of PN morphologies. 
Subject headings: stars: evolution — planetary nebulae: general — galaxies: individual (LMC) 


1. INTRODUCTION 

Despite 250 years of astronomical observations and 
decades of theoretical work, our understanding of Plan¬ 
etary Nebulae (PNe) remains rather poor. In the tra¬ 
ditional theoretical picture, the progenitors of PN e are 
low- to intermediate-mas s single stars (e.g. Abell & Gol- 
|dreich 1966; Balick 1987), which eject a large traction of 
their envelope at the end of the asymptotic giant branch 
(AGB) phase. Intense photoionization of the ejected ma¬ 
terial from the exposed stellar core then leads to the for¬ 
mation of an extended emission nebula. The widely-held 
belief that the Sun will someday become a PN is rooted 
in this picture, as is the implication that all single stars 
that go through an AGB phase will eventually form PNe. 
Yet, there is no observational evidence for this. In fact, 
the single star hypothesis cannot account for the non- 
spherical morphologies of many PNe, the observed low 
rate of PN formation per unit stellar mass, or the oc¬ 
casional appeara nce of PNe in old stellar systems like 


globular clusters (Jacoby et al. 1997; Buzzoni et al.|2006 
Moe fc De Marco|2006 ' De lVlarco|2009| ). To explain these 


these 

and other inconsistencies, new paradigms of PN forma¬ 
tion have emerg ed, invoking r otat ion and magnetic fields 
in single stars (Garcia-Segura et al. 1999), the interac¬ 
tion of AGB winds with close binary companions (Soker 


1997), and large mass transfer rates during common en¬ 
velop e episodes (Ciardullo et al. 2005 Nordhaus et al. 


2007). All of these mechanisms, or any combination of 
them, could work in principle, but without robust con¬ 
nections between PNe and their stellar progenitors, and 
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reliable measurements of mean PN lifetimes, it is diffi¬ 
cult to determine which ones are at play, and what kinds 
of PNe they produce. Unfortunately, much of what we 
know about PNe as a class comes from either case-by-case 
studies of individual objects or the analysis of small, het¬ 
erogenous samples that lack the statistical rigor needed 
to probe these issues. 

Here we present the first statistical analysis of the re¬ 
lationship between PNe and their parent stellar popula¬ 
tions in the Large Magellanic Cloud (LMC). At a well - 
determined distance of 50 kpc ( Pietrzyriski et al.] 2013), 
our closest Galactic neighbor is the ideal setting for this 
kind of study. On the one hand, the population of PNe 
in the LMC has been examined in great detail, culmi¬ 
nating in a catalog of 435 spectroscopically confirmed 
objects that extend more than 7 magnitudes down the 


[O III] A5007 PN luminosity Junction (Reid & Parker 
2010 Reid 2014, see Figure [l]). This catalog of PJNe 


is virtually free of interlopers, and is statistically com¬ 
plete to a flux limit of 10“ 14 erg cm -2 s -1 in [O III] 
A5007 (i.e., a line luminosity of 10 33 erg s -1 ) R On the 
other hand, the LMC has the best-studied stellar popu¬ 
lation in the Local Group. In particular, the stellar age 
distribution (SAD) has been mapped across the entire 
galaxy through a combination of ground-based and Hub¬ 
ble Space Telescope observations of millions of individ¬ 
ual stars, and their compa rison to theoretical isochrones 
(Harris & Zaritsky |2009 ). This SAD map consists of 
1376 cells covering the inner 64 deg 1 2 of the galaxy, with 
a spatial resolution of 12' x 12' (350 x 350 pc at 50 kpc) 
and a temporal resolution of 0.2 dex (16 bins spanning 
lookback times between 6 Myr and a Hubble time). 

This SAD map can provide unique insights into the 


4 To preserve sample homogeneity and p urity, we only use PNe 
found in field MC22 of the UKST survey (|Reid & Parker||2006|) 
whose spectra are classified as ‘true’, and discard TikelW a na h pos- 
sible’ obiects (Reid 2014), as well as outer LMC PNe (Reid~&1 
|Parker|20 13). 1 - 1 
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Figure 1 . Planetary Nebulae from||Reid & Parker||2010) (red crosses), superimposed on an HI map of the LMC (Braun 2012| grayscale). 
The area covered by the PN catalog is outlined by the orange box. 


properties of the stellar progenitors of PNe. The SAD 
is the localized version of the global star-formation his¬ 
tory, which gives the amount of stellar mass formed as a 
function of lookback time for an entire galaxy. The stars 
found today in each individual region of the LMC have 
drifted there over the years, after having been formed 
at many different times and locations throughout the 
galaxy, but the SAD of the region still provides a com¬ 
plete census, within observational uncertainties, of its 
current stellar content. As such, it must include the 
progenitors of all astronomical transients found in said 
region (including PNe), provided that the lifetimes of 
these transients are short compared to the dynamical 
drift timescale. 

2. DELAY TIME DISTRIBUTION RECOVERY METHOD 

From a statistical point of view, the relationship be¬ 
tween PNe and their stellar progenitors is encapsulated 
in the Delay Tine Distribution (DTD), defined as the 
rate of production of PNe as a function of time after a 
hypothetical brief burst of star formation. Essentially, 
the DTD is the impulse response, or Green’s function, 


of the formation rate of PNe, and as such it reveals the 
evolutionary timescales associated with PN progenitors. 
Although the DTD can be a key diagnostic for objects 
with unknown or uncertain progenitors, such as Type 
la supernovae (see Maoz et al.j 20T4]) , it is usually diffi¬ 
cult to recover from extragalactic data sets. Our work 
showcases the advantages of deriving DTDs in the Local 
Group, where we have access to resolved stellar popula¬ 
tions and highly complete object catalogs. For the case at 
hand, the current (t = to) PN formation rate, Ri(to), in 
a certain region i (i = 1... K) within the LMC is given 
by the convolution of T(£), the main-sequence turnoff 
rate of a single-age stellar population formed at time t 
(which has units of turnoff stars per year per solar mass 
formed), with the region’s SAD, rhi(t) (which has units 
of total mass of stars formed per year): 



rhi(t o — t)^(t)dt. 


(i) 


Recall that the SAD includes all stars that are now 
in region i, irrespective of where and when they were 
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formed. The number of PNe expected to be observed 
in the region, A^(to), is the product of the integrand in 
this expression and the mean lifetime of PNe, Tpw{t), 
expressed as a function of the main-sequence lifetime of 
their progenitors t: 


A. 


rto 

(to) = / ' 

J 0 


rhi(to - t)^f{t)T PN {t)dt. 


( 2 ) 


Note that Tppf(t) can equal zero if stars of a particular 
main-sequence lifetime t (i.e., formed with a particular 
mass) do not go through a PN phase. 

Given a set of observed PNe in each region NpN r i of the 
LMC, and an SAD map of the galaxy rrq(£), recovering 
the DTD (T(t)TpAr(t)) is a typical inverse problem. To 
solve it, we f ollow the procedure described in |Maoz fe; 
|Badenes| (|2010|) , re-casting the convolution in Eq. 2 as a 
discrete sum over time-binned versions of rhi(to — t ) and 
T(£)TpAr(t): 

M 

A*= y 2 rn i,j('kTp N ) j . 

3 = 1 


( 3 ) 


Here, the number of PNe expected in each spatial cell 
i of the SAD map, A i (i = 1... K) is the product of 
the discretized SAD in the cell ( rriij , the stellar mass 
formed during time interval j, with j = 1... M) and the 
discretized DTD (4/j, the number of stars that turn off 
the main sequence per unit formed stellar mass, per unit 
time during time interval j, multiplied by the mean life¬ 
time Tpnj of PNe formed from stars with main-sequence 
lifetimes in the time interval j). The values of (^Tp^)j 
are adj usted with t he Marko v Chain Monte Carlo solver 
emcee (Foreman-Mackey et al.|2013|, so that random re¬ 
alizations of a Poisson process with expectation value A i 
give, on average, the best fit to the actual number of PNe 
observed in each cell, Npnj. The value of NpN r i r anges 
between 0 and 7 in th e 732 cells in the SAD map of Har¬ 
ris & Zaritsky (2009) which overlap the area of the PN 
catalog (see Figure PH) . 

Statistical 1 and z a error intervals on each (TTptv)^ 
are calculated by applying a highest density criterion 
to the posterior probability distributions generated by 
emcee. To calculate the additional errors due to the un¬ 
certainty in the SADs, we repeat the procedure using 
the upper and lower limits to the SAD of each cell and 
temporal bin of the map, and compute the difference be¬ 
tween the new median values and those obtained from 
the best-fit solution. Statistical and SAD-related errors 
are then added in quadrature. The temporal resolution 
of the recovered DTD is determined by the size of the 
observational sample (i.e., the 435 PNe), as there is a 
fixed amount of ‘signal’ to be spread amongst the M 
temporal bins. For the LMC PNe, we have found that 
a solution with M = 6 provides the best compromise 
between temporal resolution and significance of the de¬ 
tected progenitor populations, given the statistical and 
SAD uncertainties. As a verification test, we also solved 
for (TTp 7 v)j using the adaptive grid-based search of pa¬ 
rameter space with unce rtainties based on Monte Carlo 
simulations described 
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obtained results similar to those described" 


(2010), and 
:>elow. 


The recovered DTD for the LMC PNe is presented 
in Figure [2] and Table 1. Its most prominent feature 
is the clear detection of a population of PN progenitors 
with a narrow range of main-sequence lifetimes between 
5 and 8 Gyr. These lifetimes correspond to zero-age 
main sequence masses between 1.2 and 1.0 M 0 _m_isa 


lated stella r models of LMC metallicity (Bertelli et al 
2008| |2009|). The DTD also suggests the presence of a 


second population of progenitors with main-sequence life¬ 
times between 35 and 800 Myr (masses between 8.2 and 
2.1 M 0 ), but when the SAD uncertainties are taken into 
account, the formal significance of this detection drops 
below 2cr. Nevertheless, two separate lines of evidence 
lend support to this tentative detection. First, a few well- 
obse rved LMC PNe have central core masses above ^0.7 
M 0 ( Villaver et al.| |2007), which, through the initial- to 
final-mass relation of Kalirai et al. (2008), require pro¬ 
genitors with main-sequence li fetimes below 800 M yr (in 
some cases, as low as 60 Myr; Dopita et al.||1993). Sec¬ 
ond, the spatial correlation between the distribution of 
formed stellar mass in each DTD age bin and the loca¬ 
tions of PNe in the LMC is strong enough to be seen by 
eye in both the 35-800 Myr and the 5-8 Gyr bins (see 
Figure |3|. For all other main-sequence lifetimes, the cor¬ 
relation is less clear, and we obtain only 2 a upper limits 
to (TTpat)j. 

We emphasize that these results do not depend on the 
details of the global star-formation history of the LMC 
or the SAD in each individual region, which is effectively 
a nuisance function that is disentangled from the DTD 
in our analysis. Although dynamical processes will mix 
the LMC’s stellar populati ons o n ti mescales of a few disk 

. et al.|2009), 

Tie actual die 


crossing times (^200 Myr; Bastian 


the mea¬ 
sured SAD for each region i gives the actual distribution 
of stellar ages that have produced the PNe in that specific 
location, and therefore the question of when, or where in 
the galaxy, those stars were formed is irrelevant. 

The mean lifetimes of the PNe produced in each time 
bin of the DTD can be calculated by dividing the re¬ 
covered (^Tpjsf)j for each time bin j by the rate at 
which stars from a coeval population with age t leave 
the main sequence, T(£). This rate can be derived as fol¬ 
lows. Let us assume a power-law dependence of main se¬ 
quence lifetime on initial mass, £ ms = ti(m/M 0 ) /3 , where 
t\ is the lifetime, in Gyr, of a 1 M 0 star. If the mass- 
normalized initial mass function is also a power law be¬ 
tween mass m and the minimum stellar mass, m m i n (i.e., 
dN/dm= (a + 2)(ra min /M 0 ) _ ( a+2 )(ra/M 0 ) a ), then the 
main sequence turnoff rate, in stars per Gyr per unit 
stellar mass formed, will be 


*(*) = 


a + 2 ( m 


Mr ., 


-(a+2) / t \ 


q + l 

0 


(£) 

( 4 ) 

J —2.5 and , for stars 
Gyr (Bertelli et al.| 


3. THE DTD OF PNE IN THE LMC 


From stellar evolution models, /? 
with LMC metallicities, t\ — < 

2008, 2009). If, for consistency with the SAD maps, 
we adopt a Salpeter initial mass function (a = —2.35, 
777/min = 0-1 M Q ), we obtain a turnoff rate of T(t) = 
0.020(t/Gyr) -0,46 stars per Gyr per solar mass formed. 
Dividing (4/ Tppr)j by 4/j gives the effective lifetime of 
PNe in each DTD bin. This lifetime is a mean over the 
PN population in a given age bin. If, for example, only 
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Figure 2. (a) Delay time distribution (DTD) of PNe in the LMC, showing the number of PNe produced by a single-age stellar population 
per unit formed stellar mass as a function of time since star formation. There is a clear detection in the 5 — 8 Gyr time bin, corresponding to 
main-sequence lifetimes of 1.2 — 1.0 Mq stars. A second, less significant, detection can be seen in the 35 Myr to 800 Myr bin, corresponding 
to stars with initial masses of 2.1 — 8.2 Mq. Upper limits (2cr) are shown for the other bins, (b) Mean PN lifetimes as a function of delay 
time since star formation, obtained by dividing the DTD in (a) by the rate of main sequence turnoff T(t) (Eq. 4). 


Table 1 

Delay time distribution of PNe in the LMC 


Main-sequence lifetime 

Stellar Mass a 

0 %T PN )j 

T PN,j 

[Myr] 

[M 0 ] 

[PNe / 10 6 M o ] 

[kyr] 

< 35 

> 8.2 

< 5.8 

< 14 

35 - 800 

8.2-2.1 

1 4+0-7 

H+ 6 

800 - 2000 

2.1 - 1.6 

< 0.5 

< 9 

2000 - 5000 

1.6 - 1.2 

< 0.2 

< 5 

5000 - 8000 

1.2 - 1.0 

o 

A 

1 + 
o o 

to 

27 ±6 

> 8000 

< 1.0 

< 0.1 

< 7 


a Stellar masses that corres pond to t he main-seq uenc e lifetimes listed in the first column for isolated stellar 
models of LMC metallicity ( jBert elli et al .|2008| 2 009| > 



RA (deg) RA (deg) RA (deg) 


Figure 3. PNe from|Reid Sz Parker] (|2010) (red crosses, same i n all panels), sup erimposed on the distribution of stellar mass formed in 
the six temporal bins of the DTD, taken from the SAD maps of|Harris & Zaritsky] (2009) ( gra yscale). The correlation of the PNe with the 
35-800 Myr and 5-8 Gyr stellar populations, indicated quantitatively by the DTD in Eigure[2] is discernible qualitatively by eye in the two 
middle panels. Most PNe in the LMC come from stars with main-sequence lifetimes in these two ranges. 
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half of the stars in an age bin (say, those in close binaries) 
actually produce PNe, the lifetimes of those PNe will 
be twice the value calculated in this way. This method 
yields lifetimes of 27 ± 6 kyr for the PNe produced by 
the older progenitors, and 11+f kyr for the PNe pro¬ 
duced by the younger ones (see Figure [2] and Table 1). 
These numbers are in rough agreement with estimates 
based on a combination of radiation-hydrodynamics PN 
models and observed nebular expansion velocities of local 
objects (Jacob et al. 2013). 


4. DISCUSSION AND CONCLUSIONS 

The DTD presented here allows us, for the first time, 
to examine the properties of PN progenitors in a context 
where the observational sample of PNe is highly pure 
and complete, and the entire underlying stellar popula¬ 
tion has been taken into account. This has important 
implications for our understanding of the stellar evolu¬ 
tion of low- and intermediate-mass stars. We confirm 
the theoretical expectation that some stars do not make 
a detectable contribution to the LMC PN population, 
either because they explode as core collapse supernovae 
(t < 35 Myr, M > 8.0 M 0 ) or because their post-AGB 
evolutionary timescales are longer than the timescale for 
the diss ipation of th e ejected envelope (t > 8 Gyr, M < 
1.0 M 0 ; [Herwigj 2005). Surprisingly, we find no contribu¬ 
tion to the PN population from stars with main-sequence 
lifetimes between 800 Myr and 2 Gyr, which must evolve 
through an AGB phase of some kind (Herwig 2005). If 
these stars produce PNe, they must be either extremely 
faint (with fluxes below the completeness limits of most 
PN surveys in the Milky Way), or short-lived (with mean 
lifetimes below 9 kyr), or both. 

Our DTD and derived PN lifetimes result in an in¬ 
tegrated PN formation rate of ~ 0.02 yr -1 in the sur¬ 
veyed are a, which includes ~ 80% of the stellar mass of 
the LMC ([Harris & Zaritsky||2009|). Tak ing an absolute 
magnitude of My = —18.4 for the LMC flde Va nconlenrs 
et al.|199d]) and a bolometric correction of —0.8 ( |Buzzoni 


et af 


2006), this translates into a bolometric-lummosity 

-l t -i 


specific PN formation rate of ~ 7 x 10“ 12 PNe yr 
which is lower than the typical values obtained in studies 
of Galactic PNe. However, Milky Way surveys must deal 
with heterogeneous, incomplete samples of questionable 
purity, and are affected by the systematic unc ertainties 
assoc i ated with poorly known distances ([Parker et al 


2006| |Sabin et ai.l|2014| but see |Frew et al.||2015|). Be¬ 
cause all LMC PlNe are at a well-determinecT distance, 


e are at a well-determined distance, 
and have known foreground reddenings, our analysis is 
largely free from these problems. 


are 


Among the PNe visible in the LMC today, 40^29 
generated by the younger progenitors and 25 db 8% by 
the older progenitors. Thus, roughly one third of the 
LMC PNe could come from stars in other age bins, for 
which our DTD only gives upper limits. Nevertheless, 
the presence of two distinct and disjoint progenitor pop¬ 
ulations strongly suggests the existence of two separate 
formation channels for PNe, a possi bility t hat has b een 
quali t atively discussed in the past dMoe fe De Marco 
2006; De Marco 2009 Frew & Parker 20 10|), but never 
observationally confirmed. 

If we restrict our analysis to the PNe in the bright¬ 
est third of the sample (L[ 0 m] ^ 4 x 10 34 erg s _1 ), the 


younger progenitor population is not detected. This sug¬ 
gests, somewhat counter-intuitively, that it is the older, 
less-massive progenitors that produce the brighter PNe. 
If true, this might explain a long-standing issue with the 
use of the [O III] PN luminosity function for calculating 
extragalactic distances. This method is well-calibrated 
in more than a dozen galaxie s with known Cepheid dis¬ 
tances (Ciardullo et al. 2002), yet it yields distances to 
the Virgo and Fornax clusters that are ~ 10% smaller 
than those found with o ther methods ([J acoby e t al.|1990 


McMillan et al. 1993 Freedman et al. " 20Uip . In the 


past, this offset has been attributed to the presence of an 
unidentified systema tic err or whi ch o nly affects galaxies 
beyond ~ 10 Mpc (Ferrarese et al. 2000), but a more 
physical explanation is that the brightest PNe in old 
elliptical galaxies are intrinsicall y more luminou s than 
their spiral galaxy counterparts (Ciardullo 2012). This 
is consistent with our results. 

A detailed interpretation of our DTD in the frame¬ 
work of specific PN formation mechanisms is beyond the 
scope of the present work. Nevertheless, we can still put 
together a basic picture using some additional informa¬ 
tion. The scale height of bipolar PNe in the Milky Way is 
significantly smaller than that of ro und or elliptical PNe 
(Phillips 2001 Parker et al. 2006), implying that they 
have younger progenitors (|Corradi & Schwarz||1995|). In 
the LMC, HST has imaged 68 of the 435 PlNe in our 
sample, with 49 objects showing symmetric (round or el¬ 
liptical) morph ologies a nd 19 showing distinctly bipolar 
structures (Shaw et al.|2006). The DTD of the symmetric 
subsample of PlNe shows a significant de tection in the 5- 
8 Gyr bin (as suggested by Stanghellini 2009), while the 
subsample of bipolar PNe is too small to produce any 
statistically significant detections. We have also derived 
the DTD of the 124 PNe that show enhanced N in their 
spectra (W. Re i d, priv. comm.). These ‘Type F PNe 
(|Peimbert||1978| |Peimbert fc Torres-Peimbert|)l98 3|) are 
associated wit n core masses larger than 0.64 M 0 (|Kaler 
fc Jacoby|1990 ), corresponding to main-sequence masses 


above ~2.25 M 0 . The DTD shows a marginal (1.7cr) de¬ 
tection in the 35-800 Myr bin, consistent with this min¬ 
imum core mass. Taken together, these results suggest 
a scenario wherein the 35-800 Myr progenitors produce 
fainter, more asymmetr ic PNe, possibly of binary origin 
(Nordhaus et al. 2007), while the 5-8 Gyr progenitors 
produce brighter, more symmetric PNe, perhaps via tra¬ 
ditional single star evolution. 

This basic picture is in rough agreeme nt with the prop ¬ 


erties of individual P Ne in the LMC ( Dopita et al.| 1993 
Villaver et al. 2007), and the Galactic Bulge ( |Gesicki 


et al. 2014)~but a rigorous statistical validation would 
require a complete morphological census. The main 
features of our proposed scenario are also supported 
by theoretical work. On the one hand, the PN life¬ 
times stemming from the older progenitors are consis¬ 
tent with state-of-the-art radiat ion-hydrodynamics mod¬ 
els for single-star PNe (Jacob et al.|[2013|). On the other 
hand, it has been shown that binary systems that enter 
a common envelope phase while the primary is on the 
AGB can form bipolar PNe, provided t hat th e primary 
has a minimum mass of 2 M 0 (|Soker 1998]). We leave 
a more detailed evaluation of these and other specific PN 
formation mechanisms to future work. 

We conclude by pointing out that the methods and 
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techniques presented here are applicable to a wide variety 
of astronomical objects, from variable stars (Cepheids, 
RR Lyrae, S Scuti, etc.) to interacting binaries (novae, 
cataclysmic variables, high- and low-mass X-ray bina¬ 
ries, etc.). The best available catalogs for these objects, 
and the most complete census of their progenitor stel¬ 
lar populations, are found in Local Group galaxies. The 
recovery of high-quality DTDs from these data sets can 
provide much-needed tests for specific stellar evolution 
scenarios. 
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Yam, Dan Foreman-Mackey, Dustin Lang, and Dennis 
Zaritsky. Letizia Stanghellini and Warren Reid kindly 
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of the paper, which greatly benefitted from the sugges¬ 
tions of the anonymous referee. DM acknowledges sup¬ 
port by the I-CORE Center of the Planning and Budget¬ 
ing Committee and the Israel Science Foundation. The 
Institute for Gravitation and the Cosmos is supported 
by the Eberly College of Science and the Office of the 
Senior Vice President for Research at the Pennsylvania 
State University. 
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